#######################################################
############### GOV 2001 Spring 2014 ####################
##############  Social Camouflage and the Conditional Effects of 
############# Racial Diversity on Race-Targeted Policies
# Angie M. Bautista-Chavez & Boram Lee
# Department of Government, Harvard University
#######################################################
###### For replication dataset, contact authors: boramlee01@fas.harvard.edu
#######################################################

rm(list=ls())

setwd("/Users/boramlee/Desktop")
library(foreign)
data <- read.dta("2001_Newdata.dta")

install.packages("ZeligChoice")
require(ZeligChoice)

install.packages("plm")
require(plm)

install.packages("lmtest")
require(lmtest)

install.packages("Hmisc")
require(Hmisc)

install.packages("stargazer")
require(stargazer)

install.packages("xtable")
require(xtable)

library(car)
library(rpart)





#####################################

data$age100 <-data$age/100
data$area <- rowSums(cbind(data$area, data$areacode), na.rm=TRUE)
data <- cbind(data, data$age100, data$area)

n <- nrow(data)
j <- length(unique(data$area))



########################################
############ Creating Contextual Measure. 
########################################

#################################
names(data)[names(data) == 'P002001'] <- 'population'
names(data)[names(data) == 'PCT049002'] <- 'belowpov'
names(data)[names(data) == 'HCT012001'] <- 'houseincome'
names(data)[names(data) == 'HCT027001'] <- 'housingunit'
names(data)[names(data) == 'HCT027009'] <- 'renteroccup'

############### SuPER IMPORTANT ########################
# By immobile in the code, we are measuring IMMOBILITY...
# However, we label it mobility in the paper.
data$immobile <- (data$housingunit-data$renteroccup)/data$housingunit
data$immobile <- data$immobile-mean(data$immobile)/(max(data$immobile)-min(data$immobile))
data <- cbind(data, data$immobile)


summary(data$income)
data$incomenum <- recode(data$income, "1='1'; 2='2'; 3='3'; 4='4'; 5='5'; 6='6'; 7='7'; 8='8';9='9';
                      10='10'; 11='11'; 12='12'; 13='NA'; 14='NA'")
data$incomenum <- ((as.numeric(data$incomenum) - 1)/11)
data$incomenum <- (data$incomenum-mean(data$incomenum))
hist(data$incomenum, col="magenta", border="magenta", xlim=c(-1, 1))


data$areaimmobile<-NA
data$areanew <- NA

data$bp<- data$belowpov/data$population

# Creating new area code ranging from 1 to 10.
# Replacing area codes 
for(i in 1:n){
  if(data$area[i]==212){data$areanew[i]<-1}}
for(i in 1:n){
  if(data$area[i]==718){data$areanew[i]<-2}}
for(i in 1:n){
  if(data$area[i]==914){data$areanew[i]<-3}}
for(i in 1:n){
  if(data$area[i]==845){data$areanew[i]<-4}}
for(i in 1:n){
  if(data$area[i]==516){data$areanew[i]<-5}}
for(i in 1:n){
  if(data$area[i]==631){data$areanew[i]<-6}}
for(i in 1:n){
  if(data$area[i]==518){data$areanew[i]<-7}}
for(i in 1:n){
  if(data$area[i]==607){data$areanew[i]<-8}}
for(i in 1:n){
  if(data$area[i]==315){data$areanew[i]<-9}}
for(i in 1:n){
  if(data$area[i]==716){data$areanew[i]<-10}}


# Creating weighted mean of immobility at the area level
for(i in 1:n){
  for(j in 1:10){
    if(is.na(data$areanew)==TRUE){data$areaimmobile[i]<- NA}
    else if(data$areanew[i]==j){data$areaimmobile[i] <- weighted.mean(data$immobile[data$areanew==j], data$population[data$areanew==j])}}}

# Creating weighted mean of diversity at the area level
for(i in 1:n){
  for(j in 1:10){
    if(is.na(data$areanew)==TRUE){data$areadif[i]<- NA}
    else if(data$areanew[i]==j){data$areadif[i] <- 
                                  weighted.mean(data$dif[data$areanew==j], data$population[data$areanew==j])}}}

data <- cbind(data, data$areadif, data$areaimmobile)


r <- length(unique(data$zip))
for(i in 1:n){
  for(j in 1:r){
    if(is.na(data$zip)==TRUE){data$zipnetwork[i]<- NA}
    else if(data$zip[i]==j){data$zipnetwork[i] <- mean(data$dv_1[data$zip==j], na.rm=TRUE)}}}

cbind(data$zipnetwork, data$zip, data$dv_1)


# I create mean value of dv_1 by areacode.
# It is not reported in the paper.
data$network <- NA
for(i in 1:n){
  for(j in 1:10){
    if(is.na(data$areanew)==TRUE){data$network[i]<- NA}
    else if(data$areanew[i]==j){data$network[i] <- mean(data$dv_1[data$areanew==j], na.rm=TRUE)}}}

# I create randomly sampled value of dv_1 by areacode.
# I hypothesize that a respondent's policy attitude is affected by 5 other people
# randomly selected from the area he or she is living or probably working in.
data$rnetwork <- NA
set.seed(1234)

for(i in 1:n){
  for(j in 1:10){
    if(is.na(data$areanew)==TRUE){data$network[i]<- NA}
    else if(data$areanew[i]==j){data$rnetwork[i] <- mean(sample(data$dv_1[data$areanew==j], size=10, replace=FALSE))}}}

cbind(data$areanew, data$dv_1, data$rnetwork)


################# Creating Area-Income network measure
income_bin <- NA
n <- nrow(data)

for(i in 1:n){
  if(is.na(data$incomenum[i])==TRUE){income_bin[i]<-NA}
  else if(data$incomenum[i] <= quantile(data$incomenum, .33)){income_bin[i]<-1}
  else if(quantile(data$incomenum, .33) < data$incomenum[i] && data$incomenum[i] <= quantile(data$incomenum, .66)){income_bin[i] <-2}
  else if(quantile(data$incomenum, .66) < data$incomenum[i] && data$incomenum[i] <= quantile(data$incomenum, 1)){income_bin[i]<-3}
}
data$income_bin <- income_bin




stat <- aggregate(dv_1 ~ income_bin*areanew, data=data, FUN=median)
u <- nrow(stat)

for(i in 1:n){
  for(q in 1:u){
    if(is.na(data$income_bin[i])==TRUE | is.na(data$areanew[i]==TRUE)){data$incomenet[i]<-NA}
    else if(data$areanew[i]==stat$areanew[q] && data$income_bin[i]==stat$income_bin[q])
    {data$incomenet[i] <- stat$dv_1[q]}}}

data <- data[order(data$areanew, data$libcon) , ]

cbind(data$areanew, data$income_bin, data$dv_1, data$incomenet)


################# Creating Area-wide Ideology network measure
neighbordat <- cbind(data$areanew, data$libcon, data$dv_1)

id_bin <- NA
n <- nrow(data)
quantile(data$libcon, .33)
quantile(data$libcon, .66)
quantile(data$libcon, 1)

for(i in 1:n){
  if(is.na(data$libcon[i])==TRUE){id_bin[i]<-NA}
  else if(data$libcon[i] <= quantile(data$libcon, .33)){id_bin[i]<-1}
  else if(quantile(data$libcon, .33) < data$libcon[i] && data$libcon[i] <= quantile(data$libcon, .66)){id_bin[i] <-2}
  else if(quantile(data$libcon, .66) < data$libcon[i] && data$libcon[i] <= quantile(data$libcon, 1)){id_bin[i]<-3}
}
data$id_bin <- id_bin

stat <- aggregate(dv_1 ~ id_bin*areanew, data=data, FUN=mean)

names(stat)

u <- nrow(stat)

for(i in 1:n){
  for(q in 1:u){
    if(data$areanew[i]==stat$areanew[q] && data$id_bin[i]==stat$id_bin[q])
    {data$idnet[i] <- stat$dv_1[q]}}}

data <- data[order(data$areanew, data$libcon) , ]

cbind(data$areanew, data$id_bin, data$dv_1, data$idnet)

################# Creating Area-wide Party network measure
neighbordat <- cbind(data$areanew, data$pid, data$dv_1)

pid_bin <- NA
n <- nrow(data)

for(i in 1:n){
  if(data$pid[i]==0 | data$pid[i]>0){pid_bin[i]<-1}
  else if(data$pid[i]<0){pid_bin[i] <-0}}

data$pid_bin <- pid_bin


stat <- aggregate(dv_1 ~ pid_bin*areanew, data=data, FUN=mean)
stat
names(stat)

u <- nrow(stat)

for(i in 1:n){
  for(q in 1:u){
    if(data$areanew[i]==stat$areanew[q] && data$pid_bin[i]==stat$pid_bin[q])
    {data$pidnet[i] <- stat$dv_1[q]}}}

data <- data[order(data$areanew, data$pidnet) , ]

cbind(data$areanew, data$id_bin, data$dv_1, data$pidnet)

################ Analysis starts.
############### Multi-level analysis with plm package.


############### Table 2.


require(plm)
Panel <- plm.data(data, index = c("areanew"))

#### Model 0 The Models used in Weber et al.'s paper

plm0 <- plm(formula = dv_1 ~ egalitarianism + individualism + pid + libcon + 
              gender + educ + smonitor + age100  + workbla0+ violbla0 +
              smonitor*workbla0 + smonitor*areadif + smonitor*areadif*workbla0
            + smonitor*violbla0 + smonitor*areadif + smonitor*violbla0*areadif,
            , data=Panel,
            model=c("pooling"), index=c("areanew"), family=gaussian)
summary(plm0)


# Not reported in our paper.
# Results do not change, I used pooling because variance is negative. 
# Variance of individual effects at the zipcode level is set to zero.
# Probably we do not have enough observations to estimate random effects.


# Model 1 in Table 2.
# Testing whether immobility affects people's policy attitudes.
# Main variable of interest here is a three-way interaction of dif*immobile*workbla0
plm00 <- plm(formula = dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
             + smonitor + age100  + workbla0+ incomenum + areadif+
               areaimmobile,  
             data=Panel,
             model=c("within"), index=c("areanew"), family=gaussian)

summary(plm00)
plm00<-coeftest(plm00, vcov=vcovHC(plm00, cluster="group"))
plm00


###### 

# Model 2 in Table 2.
plm01 <- plm(formula = dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
             + smonitor + age100  + workbla0+ incomenum+areadif+ 
               areaimmobile+ areadif*areaimmobile + areadif*workbla0 + workbla0*areaimmobile+
               areadif*workbla0*areaimmobile,  
             data=Panel,
             model=c("within"), index=c("areanew"), family=gaussian)

summary(plm01)
plm01_coef <-coeftest(plm01, vcov=vcovHC(plm01, cluster="group"))
plm01_coef


# Interaction of area immobility and racial diversity changes the effect of racial stereotype.
stargazer(plm00, plm01_coef, digits=2,     
          covariate.labels=c("Egalitarianism", 
                             "Individualism", 
                             "Party ID (Republican)","Ideology (Conservative)",
                             "Gender", "Education", "Self-monitor", "Age", "Lazy stereotype", "Income",
                             "Lazy x Diversity", "Lazy x Mobility", "Lazy x Diversity x Mobility"), style="apsr")


########### Figure 2 in our paper
########### Plotting Marginal Effects of Racial Stereotype
install.packages("visreg")
require(visreg)

data$area_index = factor(data$areanew)

fit <- lm(dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
          + smonitor + age100  + workbla0+ areadif+incomenum+
            areaimmobile+ areadif*areaimmobile + areadif*workbla0 + workbla0*areaimmobile+
            areadif*workbla0*areaimmobile + factor(area_index), data=data)

par(mfrow=c(2,2))

quantile(data$areadif, .1)
quantile(data$areaimmobile, .9)

# When diversity is 90% and mobility is 90%
visreg(fit, "workbla0", cond=list(areadif=.1969, areaimmobile=-.3777), xlab="Lazy stereotype", 
       ylab="Effects on Policy Opposition", type="contrast", 
       main="Diversity 90% Mobility 90%", cex.main=.8, line=list(col="blue"),
       points=list(cex=.1, pch=19))

# When diversity is high and mobility is low.
visreg(fit, "workbla0", cond=list(areadif=0.1969, 
                                  areaimmobile=.0310), xlab="Lazy stereotype", 
       type="contrast", ylab="",
       main="Diversity 90% Mobility 10%", cex.main=.8, line=list(col="blue"),
       points=list(cex=.1, pch=19))


# When diversity and mobility is high. 
visreg(fit, "workbla0", cond=list(areadif=-.06633777, 
                                  areaimmobile=-.3771907), xlab="Lazy stereotype", 
       ylab="Effects on Policy Opposition", type="contrast",
       main="Diversity 10% Mobility 90%", cex.main=.8, line=list(col="blue"),
       points=list(cex=.1, pch=19))

# When divsesity is low and mobility is low. 
visreg(fit, "workbla0", cond=list(areadif=0.1969, 
                                  areaimmobile=-.3771907), xlab="Lazy stereotype", 
       type="contrast", ylab="",
       main="Diversity 10% Mobility 10%", cex.main=.8, line=list(col="blue"),
       points=list(cex=.1, pch=19))


########### Table 3 in our paper
########### Calculating Marginal Effects
vcov=vcovHC(plm01, cluster="group")
vcov
var.beta1 <- vcov[9, 9]
var.beta4 <- vcov[11,11]
var.beta6 <- vcov[12,12]
var.beta7 <- vcov[13,13]
cov.beta1.beta4 <- vcov[9, 11]
cov.beta1.beta6 <- vcov[9, 12]
cov.beta1.beta7 <- vcov[9, 13]
cov.beta4.beta6 <- vcov[11, 12]
cov.beta4.beta7 <- vcov[11, 13]
cov.beta6.beta7 <- vcov[12, 13]

attach(data)
############ WHEN area is highly mobile and diversity is low. 
########### Racial stereotype increases opposition by 0.75. And this is significant.
am <- quantile(areaimmobile, probs=.10, na.rm=TRUE)
diverse <- quantile(areadif, probs=.10, na.rm=TRUE)
SR11 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR11, digits=3)
# -0.1

SR11.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR11.se, digits=2)
# 0.28



######################################################
# Marginal Effects of Racial Policy Opposition
# when am is 10th Percentile 
# whendiverse is 25th Percentile 
######################################################
am <- quantile(areaimmobile, probs=.1, na.rm=TRUE)
diverse <- quantile(areadif, probs=.25, na.rm=TRUE)
SR21 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR21, digits=2)
# -0.09

SR21.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR21.se, digits=2)
# 0.28

######################################################
# Marginal Effects of Racial Policy Opposition
# when am is 10th Percentile 
# when sm is 75th Percentile 
######################################################
am <- quantile(areaimmobile, probs=.1, na.rm=TRUE)
diverse <- quantile(areadif, probs=.75, na.rm=TRUE)
SR31 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR31, digits=2)
# 0.15

SR31.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR31.se, digits=2)
# 0.10

######################################################
# Marginal Effects of Racial Policy Opposition
# when am is 10th Percentile 
# when diverse is 90th Percentile 
######################################################
am <- quantile(areaimmobile, probs=.1, na.rm=TRUE)
diverse <- quantile(areadif, probs=.9, na.rm=TRUE)
SR41 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR41, digits=2)
# 0.31

SR41.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR41.se, digits=2)
# 0.04

######################################################
# Marginal Effects of Racial Policy Opposition
# when mobility is 75th Percentile 
# when diversity is 10th percentile
######################################################
am <- quantile(areaimmobile, probs=.25, na.rm=TRUE)
diverse <- quantile(areadif, probs=.10, na.rm=TRUE)
SR12 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR12, digits=2)
# 0.17

SR12.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR12.se, digits=2)
# 0.1

######################################################
# Marginal Effects of Racial Policy Opposition
# when am is 25th Percentile 
# when diverse is 25th Percentile 
######################################################

am <- quantile(areaimmobile, probs=.25, na.rm=TRUE)
diverse <- quantile(areadif, probs=.25, na.rm=TRUE)
SR22 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR22, digits=2)
# 0.17

SR22.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR22.se, digits=2)
# 0.1

######################################################
# Marginal Effects of Racial Policy Opposition
# when Mobility is 75th Percentile 
# when Diversity is 75th Percentile 
######################################################
am <- quantile(areaimmobile, probs=.25, na.rm=TRUE)
diverse <- quantile(areadif, probs=.75, na.rm=TRUE)
SR32 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR32, digits=2)
# 0.01

SR32.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR32.se, digits=2)
# 0.07

######################################################
# Marginal Effects of Racial Policy Opposition
# when Mobility is .75
# when diversity is 90th Percentile 
######################################################
am <- quantile(areaimmobile, probs=.25, na.rm=TRUE)
diverse <- quantile(areadif, probs=.9, na.rm=TRUE)
SR42 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR42, digits=2)
# -0.09

SR42.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR42.se, digits=2)
# 0.11

######################################################
# Marginal Effects of Racial Policy Opposition
# when Mobility is .25
# when Diversity is .10.
######################################################
am <- quantile(areaimmobile, probs=.75, na.rm=TRUE)
diverse <- quantile(areadif, probs=.10, na.rm=TRUE)
SR13 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR13, digits=2)
# 0.27

SR13.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR13.se, digits=2)
# 0.06

######################################################
# Marginal Effects of Racial Policy Opposition
# when Mobility is .25
# when Diversity is .25
######################################################
am <- quantile(areaimmobile, probs=.75, na.rm=TRUE)
diverse <- quantile(areadif, probs=.25, na.rm=TRUE)
SR23 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR23, digits=2)
# 0.26

SR23.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR23.se, digits=2)
# 0.06

######################################################
# Marginal Effects of Racial Policy Opposition
# when Mobility is .25
# when Diversity is .75
######################################################
am <- quantile(areaimmobile, probs=.75, na.rm=TRUE)
diverse <- quantile(areadif, probs=.75, na.rm=TRUE)
SR33 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR33, digits=2)
# -0.04

SR33.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR33.se, digits=2)
# 0.08

######################################################
# Marginal Effects of Racial Policy Opposition
# when Mobility is .25
# when Diversity is .9
######################################################
am <- quantile(areaimmobile, probs=.75, na.rm=TRUE)
diverse <- quantile(areadif, probs=.9, na.rm=TRUE)
SR43 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR43, digits=2)
# -0.24

SR43.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR43.se, digits=2)
# 0.14


######################################################
# Marginal Effects of Racial Policy Opposition
# when Mobility is .10
# when Diversity is .10
######################################################
am <- quantile(areaimmobile, probs=.9, na.rm=TRUE)
diverse <- quantile(areadif, probs=.10, na.rm=TRUE)
SR14 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR14, digits=2)
# 0.36

SR14.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR14.se, digits=2)
# 0.09

######################################################
# Marginal Effects of Racial Policy Opposition
# when Mobility is .10
# when Diversity is .25
######################################################
am <- quantile(areaimmobile, probs=.90, na.rm=TRUE)
diverse <- quantile(areadif, probs=.25, na.rm=TRUE)
SR24 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR24, digits=2)
# 0.35

SR24.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR24.se, digits=2)
# .09

######################################################
# Marginal Effects of Racial Policy Opposition
# when Mobility is .10
# when Diversity is .25
######################################################
am <- quantile(areaimmobile, probs=.9, na.rm=TRUE)
diverse <- quantile(areadif, probs=.75, na.rm=TRUE)
SR34 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR34, digits=2)
# -0.09

SR34.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR34.se, digits=2)
# 0.11

######################################################
# Marginal Effects of Racial Policy Opposition
# when Mobility is .10
# when Diversity is .90
######################################################
am <- quantile(areaimmobile, probs=.9, na.rm=TRUE)
diverse <- quantile(areadif, probs=.9, na.rm=TRUE)
SR44 <- coef(plm01)[9] + coef(plm01)[12]*am + coef(plm01)[11]*diverse + coef(plm01)[13]*am*diverse
round(SR44, digits=2)
# -0.37

SR44.se <- sqrt(var.beta1+var.beta4*(diverse^2) + var.beta6*(am^2) + var.beta7*(am^2)*(diverse^2)
                + 2*diverse*cov.beta1.beta4 + 2*am*cov.beta1.beta6 
                + 2*am*diverse*cov.beta1.beta7 + 2*am*diverse*cov.beta4.beta6
                + 2*am*(diverse^2)*cov.beta4.beta7 + 2*(am^2)*diverse*cov.beta6.beta7)
round(SR44.se, digits=2)
# 0.18

######################################################
# Make and label Table 2 (Top Half) 
######################################################
table <- matrix(c(SR11, SR21, SR31, SR41, SR11.se, SR21.se, SR31.se, SR41.se, 
                  SR12, SR22, SR32, SR42, SR12.se, SR22.se, SR32.se, SR42.se,
                  SR13, SR23, SR33, SR43, SR13.se, SR23.se, SR33.se, SR43.se,
                  SR14, SR24, SR34, SR44, SR14.se, SR24.se, SR34.se, SR44.se), ncol=8, byrow=FALSE)
colnames(table) <- c("Immobility (10th)", "", "Immobility (25th)", "(25th)", "Immobility (75th)", "(75th)", "Immobility(90th)","(90th)")
rownames(table) <- c("Diversity (10th)", "Diversity (25th)", "Diversity (75th)", "Diversity (90th)")
table <- as.table(table)
xtable(table)



######## This is Table 1 in Theory section
table <- matrix(c(1,1,1,1), ncol=2, byrow=FALSE)
colnames(table) <- c("Mobile", "Immobile")
rownames(table) <- c("Non-diverse", "Diverse")
table <- as.table(table)
xtable(table)

detach(data)


#######################################################
#######################################################

######### Hypotheses 2-1 to 2-3:
#########Table 6 in Appendix Model named "Random"
# We test whether this tells us anything. 
# Are people affected by random people in their area?
# Main variable of interest here is "rnetwork"
plm21r <- plm(formula = dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
             + smonitor + age100  + incomenum + workbla0+ areadif+areaimmobile+
               rnetwork,  
             data=Panel,
             model=c("within"), index=c("areanew"), family=gaussian)

summary(plm21r)
plm21r <-coeftest(plm21r, vcov=vcovHC(plm21r, cluster="group"))
plm21r

#### H2-2. Income network matters?
#### Variable of interest is income net.
plm22 <- plm(formula = dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
             + smonitor + age100  + incomenum + workbla0+ areadif+areaimmobile+rnetwork+
               incomenet,  
             data=Panel,
             model=c("within"), index=c("areanew"), family=gaussian)

summary(plm22)
plm22<-coeftest(plm22, vcov=vcovHC(plm22, cluster="group"))
plm22

dwtest(dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
       + smonitor + age100  + incomenum + workbla0+ areadif+areaimmobile+rnetwork+
         incomenet, data=data, alternative=c("two.sided"))

#### H2-3. Ideology network matters?
#### Variable of interest is idnet.
plm23 <- plm(formula = dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
             + smonitor + age100  + incomenum + workbla0+ areadif+areaimmobile+rnetwork+
               idnet,  
             data=Panel,
             model=c("within"), index=c("areanew"), family=gaussian)

plm23<-coeftest(plm23, vcov=vcovHC(plm23, cluster="group"))
plm23

# DW test result is 1.95 and p-value 0.42. 
dwtest(dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
       + smonitor + age100  + incomenum + workbla0+ areadif+areaimmobile+rnetwork+
         idnet, data=data, alternative=c("two.sided"))


#### H2-3. Party network matters?
#### Variable of interest is idnet.
plm24 <- plm(formula = dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
             + smonitor + age100  + incomenum+ workbla0+ areadif+areaimmobile+rnetwork+
               pidnet,  
             data=Panel,
             model=c("within"), index=c("areanew"), family=gaussian)

plm24 <-coeftest(plm24, vcov=vcovHC(plm24, cluster="group"))
plm24

# DW test result is 1.96 and p-value 0.52. 
dwtest(dv_1~egalitarianism + individualism + pid + libcon + gender + educ 
       + smonitor + age100  + incomenum+ workbla0+ areadif+areaimmobile+
         pidnet, data=data, alternative=c("two.sided"))

#### H2-4. Altogether?
#### Variable of interest is idnet.
plm25 <- plm(formula = dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
             + smonitor + age100  + incomenum + workbla0+ areadif+areaimmobile+ rnetwork +
               incomenet+ idnet + pidnet,  
             data=Panel,
             model=c("within"), index=c("areanew"), family=gaussian)

plm25 <-coeftest(plm25, vcov=vcovHC(plm25, cluster="group"))
plm25



stargazer(plm21r, plm22, plm23, plm24, plm25, digits=2,
          covariate.labels=c("Egalitarianism", 
                             "Individualism", 
                             "Party ID (Republican)","Ideology (Conservative)",
                             "Gender", "Education", "Self-monitor", "Age", "Income", "Lazy stereotype",
                          "Random Network",
                             "Economic Network", "Ideology Network", "Party Network", ""), style="apsr")




#################### SM, peer groups, and Lazy stereotype
################### This is Table 4 in our paper

# Testing how high self-monitors adapt their behavior in different networks.
# Mainly, we test whether self-monitors of certain type fit their policy attitudes
# to those of their neighbors. If so, what kind of neighbors are most important?


# Test with Randomly Selected Racial Peer groups within areas. (Model 6)
plm31r <- plm(formula = dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
             + smonitor + age100  +incomenum  + workbla0+ areadif+ areaimmobile + 
               rnetwork + rnetwork*smonitor + smonitor*workbla0 + rnetwork*workbla0 + 
               rnetwork*workbla0*smonitor,  
             data=Panel,
             model=c("within"), index=c("areanew"), family=gaussian)

plm31r<-coeftest(plm31r, vcov=vcovHC(plm31r, cluster="group"))
plm31r
# Randomly selected model is not significant although the signs are as predicted.


# Testing self-monitor's behavior in income network. (Model 7)
plm32 <- plm(formula = dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
             + smonitor + age100 + incomenum + workbla0+ areadif+ areaimmobile + 
               incomenet + incomenet*smonitor + smonitor*workbla0 + incomenet*workbla0 + incomenet*workbla0*smonitor,  
             data=Panel,
             model=c("within"), index=c("areanew"), family=gaussian)

plm32<-coeftest(plm32, vcov=vcovHC(plm32, cluster="group"))
plm32


# Testing self-monitor's behavior in ideology network. (Model 8)
plm33 <- plm(formula = dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
             + smonitor + age100  + incomenum +workbla0+ areadif+ areaimmobile + 
               idnet + idnet*smonitor + smonitor*workbla0 + idnet*workbla0 + idnet*workbla0*smonitor,  
             data=Panel,
             model=c("within"), index=c("areanew"), family=gaussian)

plm33<-coeftest(plm33, vcov=vcovHC(plm33, cluster="group"))
plm33


# Testing self-monitor's behavior in party network (Model 9).
plm34 <- plm(formula = dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
             + smonitor + age100  + incomenum +workbla0+ areadif+ areaimmobile + 
               pidnet + pidnet*smonitor + smonitor*workbla0 + pidnet*workbla0 + 
               pidnet*workbla0*smonitor,  
             data=Panel,
             model=c("within"), index=c("areanew"), family=gaussian)

plm34<-coeftest(plm34, vcov=vcovHC(plm34, cluster="group"))
plm34

###

stargazer(plm31r, plm32, plm33, plm34, digits=2,
          style="apsr")


###### Visualizing Id network. This is Figure 2 in our paper.

fit <- lm(dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
          + smonitor + age100  + incomenum +workbla0+ areadif+ areaimmobile + 
            idnet + idnet*smonitor + smonitor*workbla0 + idnet*workbla0 + 
            idnet*workbla0*smonitor + factor(area_index), 
          data=data)

par(mfrow=c(1,1))

quantile(data$libcon, .9, na.rm=TRUE)
quantile(data$idnet,.9, na.rm=TRUE)
# When idnet is the most conservative this person is surrounded by conservative people.
# and When the ideology network is highly conservatie,
# Self-monitoring scheme does not work, because even high self-monitors do not
# feel the necessity to filter out their stereotype.

################## When ID peer groups oppose the policies
visreg2d(fit, "workbla0", "smonitor", cond=list(idnet=0.611), plot.type="persp",
         xlab="Lazy stereotype", ylab="Self-monitor", zlab="Opposition to Racial policies",
         cex.lab=.7, main="", cex.main=1,zlim=c(0,1))


################## When ID peer groups support the policies
visreg2d(fit, "workbla0", "smonitor", cond=list(idnet=0.326), plot.type="persp",
         xlab="Lazy stereotype", ylab="Self-monitor", zlab="Opposition to Racial policies",
         cex.lab=.7, main="", cex.main=1, zlim=c(0,1))

################## In median-positioned ideology peer groups
visreg2d(fit, "workbla0", "smonitor", cond=list(idnet=0.4852), plot.type="persp",
         xlab="Lazy stereotype", ylab="Self-monitor", zlab="Opposition to Racial policies",
         cex.lab=.7, main="", cex.main=1, zlim=c(0,1))


################## This is for income network and self-monitors.
quantile(data$incomenum, .1)
quantile(data$incomenet, .9)
fit <- lm(dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
          + smonitor + age100  + incomenum +workbla0+ areadif+ areaimmobile + 
            incomenet + incomenet*smonitor + smonitor*workbla0 + incomenet*workbla0 + 
            incomenet*workbla0*smonitor + factor(area_index), 
          data=data)


################## When economic peer groups support the policies
visreg2d(fit, "workbla0", "smonitor", cond=list(incomenet=.2633), plot.type="persp",
         xlab="Lazy stereotype", ylab="Self-monitor", zlab="Opposition to Racial policies",
         cex.lab=.7, main="", cex.main=1, zlim=c(0,1))

################## When economic peer groups positions are median
visreg2d(fit, "workbla0", "smonitor", cond=list(incomenet=.4975), plot.type="persp",
         xlab="Lazy stereotype", ylab="Self-monitor", zlab="Opposition to Racial policies",
         cex.lab=.7, main="", cex.main=1, zlim=c(0,1))

################## When economic peer groups are opposed to policies
visreg2d(fit, "workbla0", "smonitor", cond=list(incomenet=.59125), plot.type="persp",
         xlab="Lazy stereotype", ylab="Self-monitor", zlab="Opposition to Racial policies",
         cex.lab=.7, main="", cex.main=1, zlim=c(0,1))

################# To explain further, we calculate the marginal effects of stereotype.
########### This is Table 5 in our paper.
plm33 <- plm(formula = dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
             + smonitor + age100  + incomenum +workbla0+ areadif+ areaimmobile + 
               idnet + idnet*smonitor + smonitor*workbla0 + idnet*workbla0 + idnet*workbla0*smonitor,  
             data=Panel,
             model=c("within"), index=c("areanew"), family=gaussian)

vcov=vcovHC(plm33, cluster="group")
vcov
var.beta1 <- vcov[10,10] # workbla0
var.beta4 <- vcov[13,13] # smonitor * workbla
var.beta6 <- vcov[14,14] # workbla * idnet
var.beta7 <- vcov[15,15] # workbla*idenet*smonitor
cov.beta1.beta4 <- vcov[10, 13]
cov.beta1.beta6 <- vcov[10, 14]
cov.beta1.beta7 <- vcov[10, 15]
cov.beta4.beta6 <- vcov[13, 14]
cov.beta4.beta7 <- vcov[13, 15]
cov.beta6.beta7 <- vcov[14, 15]


attach(data)

############ WHEN Selfmonitor is low (10%). and idnetwork is supportive (10%).
sm <- quantile(smonitor, probs=.10, na.rm=TRUE)
ideo <- quantile(idnet, probs=.10, na.rm=TRUE)
SR11 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR11, digits=2)
#0.48


SR11.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR11.se, digits=2)
# 0.1


############ WHEN Selfmonitor is low (10%). and idnetwork is supportive (25%).
sm <- quantile(smonitor, probs=.10, na.rm=TRUE)
ideo <- quantile(idnet, probs=.25, na.rm=TRUE)
SR21 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR21, digits=2)
#0.45


SR21.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR11.se, digits=2)
# 0.1


############ WHEN Selfmonitor is low (10%). and idnetwork is opposing. (75%).
sm <- quantile(smonitor, probs=.10, na.rm=TRUE)
ideo <- quantile(idnet, probs=.75, na.rm=TRUE)
SR31 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR31, digits=2)
#0.16


SR31.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR31.se, digits=2)
# 0.12

############ WHEN Selfmonitor is low (10%). and idnetwork is opposing. (90%).
sm <- quantile(smonitor, probs=.10, na.rm=TRUE)
ideo <- quantile(idnet, probs=.90, na.rm=TRUE)
SR41 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR41, digits=2)
#0.12


SR41.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR41.se, digits=2)
# 0.14



############ WHEN Selfmonitor is low (75%). and idnetwork is supportive (10%).
sm <- quantile(smonitor, probs=.75, na.rm=TRUE)
ideo <- quantile(idnet, probs=.10, na.rm=TRUE)
SR12 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR12, digits=2)
#0.08


SR12.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR12.se, digits=2)
# 0.08


############ WHEN Selfmonitor is low (75%). and idnetwork is supportive (25%).
sm <- quantile(smonitor, probs=.75, na.rm=TRUE)
ideo <- quantile(idnet, probs=.25, na.rm=TRUE)
SR22 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR22, digits=2)
#0.1


SR22.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR22.se, digits=2)
# 0.07


############ WHEN Selfmonitor is low (75%). and idnetwork is opposing. (75%).
sm <- quantile(smonitor, probs=.75, na.rm=TRUE)
ideo <- quantile(idnet, probs=.75, na.rm=TRUE)
SR32 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR32, digits=2)
#0.24


SR32.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR32.se, digits=2)
# 0.08

############ WHEN Selfmonitor is low (25%). and idnetwork is opposing. (90%).
sm <- quantile(smonitor, probs=.25, na.rm=TRUE)
ideo <- quantile(idnet, probs=.90, na.rm=TRUE)
coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
#0.18

sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
# 0.08


############ WHEN Selfmonitor is low (25%). and idnetwork is supportive (10%).
sm <- quantile(smonitor, probs=.25, na.rm=TRUE)
ideo <- quantile(idnet, probs=.10, na.rm=TRUE)
coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
# 0.32

sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
# 0.07



############ WHEN Selfmonitor is low (25%). and idnetwork is supportive (25%).
sm <- quantile(smonitor, probs=.25, na.rm=TRUE)
ideo <- quantile(idnet, probs=.25, na.rm=TRUE)
coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
# 0.30

sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
# 0.07



############ WHEN Selfmonitor is low (25%). and idnetwork is opposing. (75%).
sm <- quantile(smonitor, probs=.25, na.rm=TRUE)
ideo <- quantile(idnet, probs=.75, na.rm=TRUE)
coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
#0.19

sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
# 0.07


############ WHEN Selfmonitor is low (10%). and idnetwork is opposing. (90%).
sm <- quantile(smonitor, probs=.75, na.rm=TRUE)
ideo <- quantile(idnet, probs=.90, na.rm=TRUE)
SR42 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR42, digits=2)
#0.27


SR42.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR42.se, digits=2)
# 0.08




############ WHEN Selfmonitor is high (90%). and idnetwork is supportive (10%).
sm <- quantile(smonitor, probs=.90, na.rm=TRUE)
ideo <- quantile(idnet, probs=.10, na.rm=TRUE)
SR12 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR12, digits=2)
# -0.07


SR12.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR12.se, digits=2)
# 0.09


############ WHEN Selfmonitor is high (90%). and idnetwork is supportive (25%).
sm <- quantile(smonitor, probs=.90, na.rm=TRUE)
ideo <- quantile(idnet, probs=.25, na.rm=TRUE)
SR22 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR22, digits=2)
# -.04


SR22.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR22.se, digits=2)
# 0.09


############ WHEN Selfmonitor is high (90%). and idnetwork is opposing (75%).
sm <- quantile(smonitor, probs=.90, na.rm=TRUE)
ideo <- quantile(idnet, probs=.75, na.rm=TRUE)
SR32 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR32, digits=2)
#0.28


SR32.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR32.se, digits=2)
# 0.14

############ WHEN Selfmonitor is high (90%). and idnetwork is opposing (90%).
sm <- quantile(smonitor, probs=.90, na.rm=TRUE)
ideo <- quantile(idnet, probs=.90, na.rm=TRUE)
SR42 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR42, digits=2)
#0.33


SR42.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR42.se, digits=2)
# 0.15


######################################################
# Make and label Table 
######################################################
table <- matrix(c(SR11, SR11.se, SR21, SR21.se, SR31,SR31.se,  SR41,  SR41.se,
                  SR12, SR12.se, SR22, SR22.se, SR32, SR32.se, SR42, SR42.se), ncol=8, byrow=TRUE)
colnames(table) <- c("Highly Liberal (10th)", "", "Liberal (25th)", "(25th)", "Conservative (75th)", "(75th)", "Highly Conservative (90th)","(90th)")
rownames(table) <- c("SM (10th)", "SM (90th)")
table <- as.table(table)
xtable(table)

detach(data)





##### Graphics
install.packages("ggplot2")
install.packages("gcookbook")
require(ggplot2)
require(gcookbook)

# We do not use these plots in the paper though.
qplot(as.factor(data$areanew), data$dif, geom="boxplot", xlab="10 Areas", ylab="Heterogeneity")
qplot(as.factor(data$areanew), data$immobile, geom="boxplot", xlab="10 Areas", ylab="Population Immobility")

########################################
# Robustness Check:: Housing Integration
########################################

############ These are Models 3-5 in Table 2 in our paper.
# recode homeown 3 4=.
# recode integr 5 6=.
# recode locgov 5 6=.
data$homeown01
data$integr
data$locgov

data$integr <- recode(data$integr, "1='1'; 2='2'; 3='3'; 4='4'; 5='NA'; 6='NA'")

install.packages("lme4")
require(lme4)
install.packages("pglm")
require(pglm)

###### INTEGR......
Panel2 <- plm.data(data, index = c("areanew"))
pglm1 <- plm(formula = as.numeric(integr) ~ egalitarianism + individualism + pid + libcon + 
               gender + educ + smonitor + age100 + workbla0 + incomenum +areadif + 
               areaimmobile + areadif * areaimmobile + areadif * workbla0 + 
               workbla0 * areaimmobile + areadif * workbla0 * areaimmobile, 
          data=Panel2, model="within", index=("areanew"),
          , family=ordinal('logit'), R=5)

pglm1 <- coeftest(pglm1, vcov=vcovHC(pglm1, cluster="group"))
pglm1

###### LOCGOV

data$locgov <- recode(data$locgov, "1='1'; 2='2'; 3='3'; 4='4'; 5='NA'; 6='NA'")

pglm2 <- plm(formula = as.numeric(locgov) ~ egalitarianism + individualism + pid + libcon + 
               gender + educ + smonitor + age100 + workbla0 + incomenum +areadif + 
               areaimmobile + areadif * areaimmobile + areadif * workbla0 + 
               workbla0 * areaimmobile + areadif * workbla0 * areaimmobile, 
             data=Panel2, model="within", index=("areanew"),
             , family=ordinal('logit'), R=5)

pglm2 <- coeftest(pglm2, vcov=vcovHC(pglm2, cluster="group"))
pglm2

###### HOMEOWN
data$homeown <- recode(data$homeown, "1='0'; 2='1'; 3='NA'; 4='NA'")
# 0 - protection for blacks, 1- hostility

pglm3 <- plm(formula = as.numeric(homeown) ~ egalitarianism + individualism + pid + libcon + 
               gender + educ + smonitor + age100 + workbla0 + incomenum + areadif + 
               areaimmobile + areadif * areaimmobile + areadif * workbla0 + 
               workbla0 * areaimmobile + areadif * workbla0 * areaimmobile, 
             data=Panel2, model="within", index=("areanew"),
             , family=binomial('logit'), R=5)

pglm3 <- coeftest(pglm3, vcov=vcovHC(pglm3, cluster="group"))
pglm3

stargazer(plm00, plm01_coef, pglm1, pglm2, pglm3, digits=2,     
          covariate.labels=c("Egalitarianism", 
                             "Individualism", 
                             "Party ID (Republican)","Ideology (Conservative)",
                             "Gender", "Education", "Self-monitor", "Age", "Lazy stereotype", "Income",
                             "Lazy x Diversity", "Lazy x Mobility", "Lazy x Diversity x Mobility"), style="apsr")


########################################
# Robustness Check:: Economic Aid
########################################
######### These are not reported in our paper.

###### HELPBLC Government general aid......
data$helpblc
data$helpblc <- recode(data$helpblc, "1='0'; 2='1'; 3='NA'; 4='NA'")
# 1 is favorable, 2 is hostility

pglm4 <- plm(formula = as.numeric(helpblc) ~ egalitarianism + individualism + pid + libcon + 
               gender + educ + smonitor + age100 + workbla0 + areadif + 
               areaimmobile + areadif * areaimmobile + areadif * workbla0 + 
               workbla0 * areaimmobile + areadif * workbla0 * areaimmobile, 
             data=Panel2, model="within", index=("areanew"),
             , family=binomial('logit'), R=5)


coeftest(pglm4, vcov=vcovHC(pglm4, cluster="group"))
pglm4

###### SPENDED
data$spended
data$spended <- recode(data$spended, "1='1'; 2='2'; 3='3'; 4='4'; 5='NA'; 6='NA'")
# 4 is most hostile

pglm5 <- plm(formula = as.numeric(spended) ~ egalitarianism + individualism + pid + libcon + 
               gender + educ + smonitor + age100 + workbla0 + areadif + 
               areaimmobile + areadif * areaimmobile + areadif * workbla0 + 
               workbla0 * areaimmobile + areadif * workbla0 * areaimmobile, 
             data=Panel2, model="within", index=("areanew"),
             , family=ordinal('logit'), R=5)


pglm5 <- coeftest(pglm5, vcov=vcovHC(pglm5, cluster="group"))
pglm5

####### Affirmative actions
###### AFFACT1
data$affact1
data$affact1 <- recode(data$affact1, "1='1'; 2='2'; 3='NA'; 4='NA'")
# 2 is most hostile

pglm6 <- plm(formula = as.numeric(affact1) ~ egalitarianism + individualism + pid + libcon + 
               gender + educ + smonitor + age100 + workbla0 + areadif + 
               areaimmobile + areadif * areaimmobile + areadif * workbla0 + 
               workbla0 * areaimmobile + areadif * workbla0 * areaimmobile, 
             data=Panel2, model="within", index=("areanew"),
             , family=binomial('logit'), R=5)
pglm6
summary(pglm6)
coeftest(pglm6, vcov=vcovHC(pglm6, cluster="group"))

###### AFFACT2
data$affact2
data$affact2 <- recode(data$affact2, "1='1'; 2='2'; 3='NA'; 4='NA'")
# 2 is most hostile

pglm7 <- plm(formula = as.numeric(affact2) ~ egalitarianism + individualism + pid + libcon + 
               gender + educ + smonitor + age100 + workbla0 + areadif + 
               areaimmobile + areadif * areaimmobile + areadif * workbla0 + 
               workbla0 * areaimmobile + areadif * workbla0 * areaimmobile, 
             data=Panel2, model="within", index=("areanew"),
             , family=binomial('logit'), R=5)
pglm7
summary(pglm7)
coeftest(pglm7, vcov=vcovHC(pglm7, cluster="group"))



############ Bottom of Table 5. Simple Slopes of Marginal Effects of negative stereotype on Policy attitudes
############ As Economic peer groups' opinions change and Self-monitoring orientations change.

plm32 <- plm(formula = dv_1 ~ egalitarianism + individualism + pid + libcon + gender + educ 
             + smonitor + age100 + incomenum + workbla0+ areadif+ areaimmobile + 
               incomenet + incomenet*smonitor + smonitor*workbla0 + incomenet*workbla0 + incomenet*workbla0*smonitor,  
             data=Panel,
             model=c("within"), index=c("areanew"), family=gaussian)

vcov=vcovHC(plm32, cluster="group")
vcov
var.beta1 <- vcov[10,10] # workbla0
var.beta4 <- vcov[13,13] # smonitor * workbla
var.beta6 <- vcov[14,14] # workbla * incomenet
var.beta7 <- vcov[15,15] # workbla*idenet*smonitor
cov.beta1.beta4 <- vcov[10, 13]
cov.beta1.beta6 <- vcov[10, 14]
cov.beta1.beta7 <- vcov[10, 15]
cov.beta4.beta6 <- vcov[13, 14]
cov.beta4.beta7 <- vcov[13, 15]
cov.beta6.beta7 <- vcov[14, 15]


attach(data)

############ WHEN Selfmonitor is low (10%). and incomenetwork is supportive (10%).
sm <- quantile(smonitor, probs=.10, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.10, na.rm=TRUE)
coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
#0.55

sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
     + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
     + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
     + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
# 0.17


############ WHEN Selfmonitor is low (10%). and incomenetwork is supportive (25%).
sm <- quantile(smonitor, probs=.10, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.25, na.rm=TRUE)
coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo

#0.39

sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
     + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
     + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
     + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)

# 0.09


############ WHEN Selfmonitor is low (10%). and incomenetwork is opposing. (75%).
sm <- quantile(smonitor, probs=.10, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.75, na.rm=TRUE)
coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
# 0.22

sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
     + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
     + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
     + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
# 0.10


############ WHEN Selfmonitor is low (10%). and incomenetwork is opposing. (90%).
sm <- quantile(smonitor, probs=.10, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.90, na.rm=TRUE)
coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
# 0.14



SR41.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR41.se, digits=2)
# 0.13



############ WHEN Selfmonitor is low (75%). and incomenetwork is supportive (10%).
sm <- quantile(smonitor, probs=.75, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.10, na.rm=TRUE)
SR12 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR12, digits=2)
#0.04


SR12.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR12.se, digits=2)
# 0.14


############ WHEN Selfmonitor is low (75%). and incomenetwork is supportive (25%).
sm <- quantile(smonitor, probs=.75, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.25, na.rm=TRUE)
SR22 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR22, digits=2)
#0.13


SR22.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR22.se, digits=2)
# 0.06


############ WHEN Selfmonitor is low (75%). and incomenetwork is opposing. (75%).
sm <- quantile(smonitor, probs=.75, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.75, na.rm=TRUE)
SR32 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR32, digits=2)
#0.21


SR32.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR32.se, digits=2)
# 0.09

############ WHEN Selfmonitor is low (25%). and incomenetwork is opposing. (90%).
sm <- quantile(smonitor, probs=.25, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.90, na.rm=TRUE)
coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
#0.19

sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
     + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
     + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
     + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)



############ WHEN Selfmonitor is low (25%). and incomenetwork is supportive (10%).
sm <- quantile(smonitor, probs=.25, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.10, na.rm=TRUE)
coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
# 0.35

sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
     + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
     + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
     + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
# 0.15



############ WHEN Selfmonitor is low (25%). and incomenetwork is supportive (25%).
sm <- quantile(smonitor, probs=.25, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.25, na.rm=TRUE)
coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
# 0.28

sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
     + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
     + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
     + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
# 0.07



############ WHEN Selfmonitor is low (25%). and incomenetwork is opposing. (75%).
sm <- quantile(smonitor, probs=.25, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.75, na.rm=TRUE)
coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
#0.22

sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
     + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
     + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
     + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
# 0.07


############ WHEN Selfmonitor is low (10%). and incomenetwork is opposing. (90%).
sm <- quantile(smonitor, probs=.75, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.90, na.rm=TRUE)
SR42 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR42, digits=2)
#0.25


SR42.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR42.se, digits=2)
# 0.13




############ WHEN Selfmonitor is high (90%). and incomenetwork is supportive (10%).
sm <- quantile(smonitor, probs=.90, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.10, na.rm=TRUE)
SR12 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR12, digits=2)
# -0.16


SR12.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR12.se, digits=2)
# 0.16


############ WHEN Selfmonitor is high (90%). and incomenetwork is supportive (25%).
sm <- quantile(smonitor, probs=.90, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.25, na.rm=TRUE)
SR22 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR22, digits=2)
# 0.02


SR22.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR22.se, digits=2)
# 0.09


############ WHEN Selfmonitor is high (90%). and incomenetwork is opposing (75%).
sm <- quantile(smonitor, probs=.90, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.75, na.rm=TRUE)
SR32 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR32, digits=2)
#0.21


SR32.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR32.se, digits=2)
# 0.13


############ WHEN Selfmonitor is high (90%). and incomenetwork is opposing (90%).
sm <- quantile(smonitor, probs=.90, na.rm=TRUE)
ideo <- quantile(incomenet, probs=.90, na.rm=TRUE)
SR42 <- coef(plm33)[10] + coef(plm33)[13]*sm + coef(plm33)[14]*ideo + coef(plm33)[15]*sm*ideo
round(SR42, digits=2)
#0.3


SR42.se <- sqrt(var.beta1+var.beta4*(sm^2) + var.beta6*(ideo^2) + var.beta7*(sm^2)*(ideo^2)
                + 2*sm*cov.beta1.beta4 + 2*ideo*cov.beta1.beta6 
                + 2*sm*ideo*cov.beta1.beta7 + 2*sm*ideo*cov.beta4.beta6
                + 2*ideo*(sm^2)*cov.beta4.beta7 + 2*(ideo^2)*sm*cov.beta6.beta7)
round(SR42.se, digits=2)
# 0.18



####################### FIN!!!!!

